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Abstract 

Hessian recovery has been commonly used in mesh adaptation for obtaining the required 
magnitude and direction information of the solution error. Unfortunately, a recovered Hessian 
from a linear finite element approximation is nonconvergent in general as the mesh is refined. It 
has been observed numerically that adaptive meshes based on such a nonconvergent recovered 
Hessian can nevertheless lead to an optimal error in the finite element approximation. This also 
explains why Hessian recovery is still widely used despite its nonconvergence. In this paper we 
develop an error bound for the linear finite element solution of a general boundary value problem 
under a mild assumption on the closeness of the recovered Hessian to the exact one. Numerical 
results show that this closeness assumption is satisfied by the recovered Hessian obtained with 
commonly used Hessian recovery methods. Moreover, it is shown that the finite element error 
changes gradually with the closeness of the recovered Hessian. This provides an explanation on 
how a nonconvergent recovered Hessian works in mesh adaptation. 
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1 Introduction 



Gradient and Hessian recovery has been commonly used in mesh adaptation for the numerical 
solution of partial differential equations (PDEs); e.g., see (3j |4j [l2j [TBJ [2lJ [23j [24] . The use typically 
involves the approximation of solution derivatives based on a computed solution defined on the 
current mesh (recovery), the generation of a new mesh using the recovered derivatives, and the 
solution of the physical PDE on the new mesh. These steps are often repeated several times until 
a suitable mesh and a numerical solution defined thereon are obtained. As the mesh is refined, 
a sequence of adaptive meshes, derivative approximations, and numerical solutions results. A 
theoretical and also practical question is whether this sequence of numerical solutions converges to 
the exact solution. Naturally, this question is linked to the convergence of the recovered derivatives 
used to generate the meshes. It is known that recovered gradient through the least squares fitting 
or polynomial preserving techniques |2l] is convergent for uniform or quasi-uniform meshes 



or 



and superconvergent for mildly structured meshes [20] as well for a type of adaptive mesh 19 
the Hessian, it has been observed that, unfortunately, a convergent recovery cannot be obtained from 
linear finite element approximations for general nonuniform meshes [2 14 , although Hessian recovery 



is known to converge when the numerical solution exhibits superconvergence or supercloseness for 



some special meshes 5 6, 15 



Surprisingly, numerical experiments also show that the numerical solution obtained with an 
adaptive mesh generated using a nonconvergent recovered Hessian is often not only convergent but 
also has an error comparable to that obtained with the exact analytical Hessian. To demonstrate 
this, we consider a Dirichlet boundary value problem (BVP) for the Poisson equation 

f-A« = /, inO = (0,1) x (0,1) 
I u = g, on cW7 

where / and g are chosen such that the exact solution of the BVP is given by 

u(x,y) = x 2 + 25y 2 . (2) 

Two Hessian recovery methods, QLS (quadratic least squares fitting) and WF (weak formulation) 
are used. (See Sect. [4] for the description of these and other Hessian recovery techniques). Figure [T] 
shows the error in recovered Hessian and the linear finite element solution with exact and recovered 
Hessian. One can see that the finite element error is convergent and almost undistinguishable for 



the exact and approximate Hessian (Fig. la) whereas the error of the Hessian recovery remains 0(1) 
(Fig. lb). Obviously, this indicates that a convergent recovered Hessian is not necessary for the 
purpose of mesh adaptation. Of course, an irrelevant recovered Hessian does not serve the purpose 
either. 

How accurate should a recovered Hessian be for the purpose of mesh adaptation? This issue has 
been studied by "Adaptive generation of quasi-optimal tetrahedral meshes" [I] and "An adaptive 
algorithm for quasioptimal mesh generation" [18]. Particularly, they show that a mesh based on an 
approximation (denoted by R) of the Hessian (denoted by H) is quasi-optimal if there exist small 
positive numbers e, 5 such that 

max\\H(x) --ff^lU ^ SX min (R(xi)), (3) 

\\R(xi) - fluj^ < eA min ( J R( : r,)) (4) 

hold for any mesh vertex x% and its patch u)i, where H Ui is the Hessian at a point in u>i where 
|det H(x)\ attains its maximum and A m i n (-) denotes the minimum eigenvalue of a matrix. Notice 
that Q does not require R to converge to H as the mesh is refined. Instead, it requires the 
eigenvalues of to be around one (cf. Sect. 0). Unfortunately, it is still too restrictive to 
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Figure 1: Error of the finite element solution and recovered Hessian as function of the number of 
mesh elements. 



be satisfied by most examples we tested; see Sect. |5j Thus, the work [T 18 does not give a full 
explanation why a nonconvergent recovered Hessian works in mesh adaptation. 

The objective of the paper is to present a study on this issue. To be specific, we consider a BVP 
and its linear finite element solution with adaptive meshes generated from a recovered Hessian. 



We adopt the M- uniform mesh approach |10 12 to view any adaptive mesh as a uniform one in 



some metric. An advantage of the approach is that the relation between the recovered Hessian 
and a mesh generated using it can be fully characterized through the so-called alignment and 
equidistribution conditions (e.g., see ([7]) and Q in Sect. [2]). This characterization plays a crucial 
role in the development of a bound for the H semi-norm of the finite element error. The bound 
converges at a first order rate in terms of the average element diameter, N~d l where N is the number 
of elements and d is the dimension of the physical domain. Moreover, the bound is valid under 
a condition on the closeness of the recovered Hessian to the exact one; see (16) or (29) and (30). 
This closeness condition is much weaker than Q. Roughly speaking, Q requires the eigenvalues 
of R~ l H to be around one whereas the new condition only requires them to be bounded below 
from zero and from above. Numerical results in Sect. [5] show that the new closeness condition is 
satisfied by all examples for four commonly used Hessian recovery techniques considered in this 
paper whereas Q is satisfied only by some examples. Furthermore, the error bound is linearly 
proportional to the ratio of the maximum (over the physical domain) of the largest eigenvalues of 
R^H to the minimum of the smallest eigenvalues. Since the ratio is a measure of the closeness of 
the recovered Hessian to the exact one, the linear dependence implies that the finite element error 
changes gradually with the closeness of the recovered Hessian. Hence, the error for the linear finite 
element approximation of the BVP is convergent for the considered Hessian recovery techniques and 
insensitive to the closeness of the recovered Hessian to the exact one. This provides an explanation 
how a nonconvergent recovered Hessian works for mesh adaptation. 

An outline of the paper is as follows. Convergence analysis of the linear finite element approxima- 
tion is given in Sects. [2] and [3] for the cases with positive definite and general Hessian, respectively. 
A brief description of four common Hessian recovery techniques is given in Sect. [4] followed by 
numerical examples in Sect. [5} Finally, Sect. [6] contains conclusions and further comments. 
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2 Convergence of linear finite element approximation for positive 
definite Hessian 



We consider the BVP 

in 0, 




(5) 

on oil 

where is a polygonal or polyhedral domain of W 1 (d> 1), £ is an elliptic second-order differential 
operator, and / and g are given functions. We are concerned with the adaptive mesh solution of 
this BVP using the conventional linear finite element method. Denote a family of simplicial meshes 
for f2 by {7~h} and the corresponding reference element by K. For each mesh Th, we denote the 
corresponding finite element solution by Uh- Cea's lemma implies that the finite element error is 
bounded by the interpolation error, i.e., 

\u - u h \ H i {n) < C\u - Uu\ H i^ n) , (6) 

where C is a constant independent of u and Th and II is the interpolation operator associated with 
the linear finite element space defined on Th- Note that ^ is valid for any mesh. 

In this paper we consider adaptive meshes generated based on a recovered Hessian R and use the 
M-uniform mesh approach with which any adaptive mesh is viewed as a uniform one in some metric 
M (defined in terms of R in our current situation) . It is known |10[ |12| that such an M- uniform 
mesh satisfies the equidistribution and alignment conditions, 



\K\det(M K )^ = ^ \K\det(M R )K VKe% (7) 

K&T h 

^tr((F' K ) T M K F' K )=det((F' K ) T M K F' K )K VK G T h (8) 

where N is the number of mesh elements, Mk is an average of M over K, Fx '■ K — > K is the 
affine mapping from the reference element K to a mesh element K, F' K is the Jacobian matrix of 
Fx (which is a constant matrix on K), and det(-) and tr(-) denote the determinant and trace of a 
matrix, respectively. In practice, it is more realistic to generate less restrictive quasi- M-uniform 
meshes that satisfy 

AT I TS\ T\/T,,\ h 

< C eq , yK e T h (9) 



N\K 


det(M K ) 2 




k 


det(M R )^ 



\tr[{F' K ) T M K F> K 
\K\^det{M K )d 



<C ali , VKe% (10) 



where C eq , C a u > 1 are some constants independent of K, N, and Th- Notice that conditions @-([8]) 
and (|9|-(10) are equivalent if C eq = C a u = 1. Numerical experiments (cf. JXTJ and Sect. [5]) show 
that quasi- M-uniform meshes with relatively small C eq and C a u can be generated in practice (cf. 
Figs. [2] and 4a). For this reason, we use quasi- M-uniform meshes in our analysis and numerical 



experiments. 

We want to investigate the effects of the closeness of the recovered Hessian R to the exact Hessian 
H on the linear finite element error with quasi- M-uniform meshes. We consider in this section a 
special case where the Hessian of the solution is uniformly positive definite on 0; i.e., there exists a 
constant 7 > such that 

H(x) > 7I, \/xen (11) 



4 



where the greater-than-or-equal sign means that the difference between the ieft-hand side and 
right-hand side terms is positive semidefinite. We also assume that the recovered Hessian R is 
uniformly positive definite on Q. This assumption is not essential and will be dropped for the 
general situation which we shall discuss in the next section. 

Recall from ^ that the finite element error is bounded by the H 1 semi-norm of the interpolation 



error of the exact solution. From 11 a metric tensor for the H 1 semi-norm can be defined as 



M K = det(R K y^2\\R K \\f 2 R K , VKeT h 



(12) 



where Rk is an average of R over K. For this metric tensor, mesh conditions ([9]) and (10) become 

(13) 



N\K\det(R K )^ \\RkM 



d 
d+2 



— <C eq , VKeTh 



E^|^|det(^)^||^|| 2 d+2 
^tv((F' K ) T R K F' K 
| if | a det(R K )* 



<C ali , VK£T h 



(14) 



Moreover, alignment condition (|14|) implies the inverse alignment condition 



\tr {{F' K )- T R^jF^)- 1 
\K\~tl det(R K )~^ 



< {dC ali ) d - ] 



(15) 



Theorem 2.1 (Positive definite Hessian). Assume that H(x) and the recovered Hessian R are 
uniformly positive definite on £1. Assume also that R satisfies 

Cr-I < R K l H{x) < Cb+I, Vx eK, MK € T h (16) 

for some positive constants Cr+ and Cr- independent of the mesh. Then, for any quasi-M -uniform 
mesh satisfying (13) and (14) associated with the metric tensor (12), the linear finite element error 



for BVP ^ is bounded by 



\ u ~ u h\H 1 (n) 



< 



Cr+ 
AT3 Cr. 



d+2 

f~< 2 (~i 2d 



det(H) • H 



Ld+2 (n) 



(17) 



Proof. From |12[ Theorem 5.1.5], the interpolation error on K is bounded by 



\u - Ii K u\ m ^ K) < C (F' h 



K 



(F' K f\H(x)\F' K 



dx 



(18) 



where |i?(a;)| = y H(x) and is the restriction of interpolation operator n on K. Notice that 
|if(a;)| = H{x) for the current situation where H{x) is symmetric and positive definite. Moreover, 
we have 



{F' K )~ l n = {F' K )~ T (F' K )~ l 



< 



R K {F' h 



K) 



\Rk\ 



Combining the above results with ([6| and using (13), ( |14[ ) and (15), we have 



dx 



\u-u h \ 2 Hl{Q) <C^\\(F' K )- 1 f [ \{F' K ) T \H{x)\F', 

K JK 

<C^kF' K )- T RJ c \F , K )-\\\R K \\ 2 - f Uf , k ) T R k F'JI\\k£h(x) 



dx 
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< CC tu E \ K V l det(R K )-^ \\Rk\\ 2 ■ C 2 aH \K\Uet(R K y ■ \K\ C 2 R+ 

K 

< CC d + l C 2 R+ • £ \Kf? det(^)3 \\R K \\ 2 



A 



d+2 



= cc^ci^imd^RK)^ \\R K \\t 

K V 

< CC^Cl + Y. (c eq ^Y,\K\det(R k )^ \\Rs\\F 

k 

\J2\K\det(R K )^ \\R K \\t^ 

d+2 2 I f 1 

CC d ^C\+C4 N- d (J2 det(R K )^\\R K 

\ K JK 



K 

d+2 2 



d+2 



d + 2 
d 



d \ d 

d 2 +2 dx) 



Assumption (PL6|) implies 



det(R K ) < det(H(x)) H-\x)R K < C]£ det(H(x)), 



\\ r k\\ 2 
Thus, we have 

I I 2 



H(x)H- l (x)R K < ||fT(aj)|| a H- L (x)R K 



< Cr- \\H{x) 



(19) 
(20) 



d+2 



det(H(x))\\H(x)\\ 2 ) dxj 

d+2 



d 
d+2 



det(#(a;)) • H(x) 



dx 



which gives (17). 



□ 



This theorem shows how a nonconvergent recovered Hessian works in mesh adaptation. Indeed, 
from (|17|) we can see that the error bound is linearly proportional to the ratio Cr+/Cr-. Since 



Cr+ I Cr- is a measure for the closeness of R to H, (17) implies that the finite element error changes 



gradually with the closeness of the recovered Hessian. 

If it! is a good (but not necessarily convergent) approximation to H, we have Cr + /Cr- = 0(1). 
Then, the error bound has a solution-dependent factor 



det(H) • H 



(21) 



Ld+2 (fi) 



On the other hand, if R is not a good approximation to H, the error will have a larger solution- 
dependent factor. For example, consider the case of R being the identity matrix which leads to the 



uniform mesh refinement. Notice that the condition (16) is mathematically equivalent to 



R^H(x) < C R+ , H-\x)R K 



1 

2 - Cr~- 



Vx 6 K, MK € %■ 



(22) 



For R = I these conditions are satisfied with 



C R+ = max ||iT(a;)||2 = max X mauX (H(x)), 

x6S2 a;Gf2 
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C R - = max ||.ff (sc)|| 2 = min X min (H(x)), 

where A max (i?(a;)) and AmiTi (H(x)) denote the maximum and minimum eigenvalues of H(x), 
respectively. Thus, for the case R = I the error bound contains a solution-dependent factor 



maxggn X ma , x (H(x)) 
m.in xen X min (H(x)) 



det(H) ■ H 



(ft) 



which is obviously larger than (21). 

Next, we study the relation between (16) and Q. Equivalence of matrix norms implies that Q 
is equivalent to 

\\R K -H{x)\\ 2 <eX min {R K ), VxeK, VKeT h (23) 



where e is a small positive constant. The above inequality yields 



R£H(x) - I < e, Vcc G K, \/K G Th 



(24) 



which in turn implies (16) with Cr+ = 1 + e and Cr- = 1 — e if e < 1. Obviously, conditions (24) 
and therefore Q require the eigenvalues of Rj c 1 H(x) to stay closely around one. On the other hand, 
condition ( |16[ ) only requires the eigenvalues of Rj c 1 H(x) to be bounded from above and below from 
zero, which is much weaker than Q. 

If R converges to H(x), both Ml) and (16) can be satisfied. However, if R does not converge to 
H(x), as is the case for most adaptive computation, the situation is different. As we shall see in 
Sect. [5j condition (16) is satisfied for all of the examples tested whereas condition Q is satisfied 
only for some examples. 

Furthermore, we would like to point out that it is unclear if (17) is optimal. However, the bound 
seems to be the best we can get. For example, if we choose the monitor function to be 

M K = det{R K )~^R K , \/K G % (25) 



which is optimal for the L 2 norm (cf. 11 ), the error bound becomes 
\ u - n /il#i(n) ^ 



C Cf/4- ^±1 i±± 



Nd Cr. 



J ali ^ e <l 



det(H) 



2 

d+4 
2d 

Ld+4. (O) 



d+4 



det(iT) • H 



(26) 



This bound has a larger solution-dependent factor than (17) since Holder's inequality yealds 



\jdet{H) ■ H 


1 
2 

< 


tfdet(H) 


2 

d+4, 
2d 




Ld+2 (n) 




Ld.+4(U) 



d+4 



det(H) • H 



Finally, it is worth mentioning that when the metric tensor (25) is used, the I? norm of the 
piecewise linear interpolation error is bounded by 



1 tt n C Cr+- 

\u-nu\\ L 2 {n) < —2 ■ - r - 



d+4 

CaliCeq 



det(H) 



2d 

Ld+4 (fi) 



(27) 



This bound is optimal in terms of convergence order and solution-dependent factor; e.g., see (7, 13 



7 



3 Convergence of the linear finite element approximation for a general 
Hessian 



In this section we consider the general situation where H(x) is symmetric but not necessarily 
positive definite. In this case, it is unrealistic to require the recovered Hessian R to be positive 
definite. Thus, we cannot use R directly to define the metric tensor which is required to be positive 
definite. A commonly used strategy is to replace R by \R\ = \^R? since \R\ retains the eigensystem 
of R. However, \R\ can become singular locally. To avoid this difficulty, we regularize \R\ with a 
regularization parameter a\ t > (to be determined). 



From (12), we define the regularized metric tensor as 
M K = det(a h I + \R K 



d+2 



\a h I + \R K \\\i +2 {a h I + \R K \), VK e T h 
and obtain the following theorem with a similar proof to that of Theorem |2.1[ 



(28) 



Theorem 3.1 (General Hessian). For a given positive parameter > 0, we assume that the 
recovered Hessian R satisfies 

(a h I + \Rk\V 1 \H(x)\ < C R+ I, Vx e K, E T h (29) 

(a h I+\R K \)- l (a h I + \H(x)\)>C R -I, \/x E K, VK (30) 

for some positive constants Cr+ and Cr- independent of the mesh. Then, for any quasi-M -uniform 
mesh satisfying ([9| and (10) associated with metric tensor (28), the linear finite element error for 
BVP ([5]) is bounded by 



\ u ~ u h\m{n) 



< 



C Cr + 
Cr- 



d+l 
1 2 

'all 



d+2 
< 2< 
'eq 



det(a h I+ \H\) ■ {a h I + \H\ 



1 

\ • (31) 



From (29), (30), and (31), we see that the greater ah is, the easier the recovered Hessian satisfies 



(29) and (30); however, the error bound increases as well. For example, consider the extreme case 
of ah —> 00. In this case, (29) and (30) can be satisfied with Cr+ = Cr- = 1 for any R. In the 



same time, the metric tensor defined in (28) has an asymptotic behavior Mk —> cn h + I and any 
corresponding M-uniform mesh is a uniform mesh. Obviously, the right-hand side of (31) is large 
for this case. Another extreme case is the limit ah — > where (31 ) reduces to (17) if both R and 
H{x) are positive definite. 



We now consider the choice of ah- It is easy to show that the equation 



tj&et(al +\H\) ■ (al + \H\) 



Ld+2 (n) 

has a unique solution a > when ^det(|i?|) • H 
element error is bounded by 

C Cr + 



Ld+2 (n) 



det(\H\)-H 

"Ld+2 (fi) 

> 0. If we choose ah 



\ u u h\H 1 (Q) 



< 



iVa Cr. 



d+l d+2 
f~< 2 (~1 2d 

u ali ° e 



d Jdet{\H\)-H 



Ld+2 (CI) 



(32) 

q, the finite 
(33) 



which is essentially the same as (17). Note that (32) is impractical since it requires the prior 
knowledge of H{x). In practice it can be replaced by 



J2\K\ det (a h I + \R K \) d + 2 \\a h I+\R K \ 

K 



d 

I d+2 
12 



K 



1 

d+2 



\RkM +2 



(34) 



Numerical results show that ah is close to a (see Fig. 4d). 
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4 A selection of commonly used Hessian recovery methods 



In this section we give a brief description of four commonly used Hessian recovery algorithms for 
two-dimensional mesh adaptation. The interested reader is referred to [14 17 for a more detailed 
description of these Hessian recovery techniques. 

Recall that the goal of Hessian recovery in the current context is to find approximations of the 
Hessian at mesh nodes using the linear finite element solution Uh- The approximation of the Hessian 
on an element is calculated as the average of the nodal approximations of the Hessian at the vertices 
of the element. 



4.1 QLS: quadratic least squares fitting to nodal values 

This method involves the fitting of a quadratic polynomial to nodal values of Uh at a selection of 
neighboring nodes in the least square sense and subsequent differentiation. The original purpose 
of the QLS was the gradient recovery (e.g., see "A new finite element gradient recovery method: 
superconvergence property" (2l]). However, it is easily adopted for the Hessian recovery by simply 
differentiating the fitting polynomial twice. 

More specifically, for a given node (say xq) at least five neighboring nodes are selected. A 
quadratic polynomial (denoted by p) is found by least squares fitting to the values of at the 
selected nodes. The linear system associated with the least squares problem usually has full rank 
and a unique solution. If it does not, additional nodes from the neighborhood of xq are added to 
the selection until the system has full rank. An approximation to the Hessian of the solution u at 
Xq is defined as the Hessian of p, viz., 

RQ LS (x ) = H(p)(x ). 



4.2 DLF: double linear least squares fitting 

The DLF method computes the Hessian by using linear least squares fitting twice. First, the least 
squares fitting of the nodal values of Uh in a neighbourhood of xq is employed to find a linear fitting 
polynomial p. The recovered gradient of function u at xq is defined as the gradient of p at xq, i.e., 

V% LF u(x ) = Vp(ajo). 

Second-order derivatives are then obtained by subsequent application of this linear fitting to the 
calculated first order derivatives. Mixed derivatives are averaged in order to obtain a symmetric 
recovered Hessian. 



4.3 LLS: linear least squares fitting to first-order derivatives 

The LLS method is similar to DLF except that the first-order derivatives at nodes are calculated in 
a different way. In this method, the first-order derivatives are first calculated at element centers 
and then at nodes by linear least squares fitting to their values at element centers. 



4.4 WF: weak formulation 

This approach recovers the Hessian by means of a variational formulation ||8|. More specifically, let 
4>o be a canonical piecewise linear basis function at node xq. Then the nodal approximation u xx ^ 
to the second-order derivative u xx at Xi is defined through 

u xx ,h{xo) [ <po{x) dx = - [ dx. 
Jn Jn ox dx 
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The same approach is used to compute u X y,h an d u yyy h- Since 0o are piecewise linear and vanish 
outside the patch associated with xq, the involved integrals can be computed efficiently with 
appropriate quadrature formulas over a single patch. 



5 Numerical examples 



In this section we present some numerical results obtained for two examples to verify the analysis 



given in the previous sections. Special attention will be paid to mesh conditions (|9|) and (10) and 



closeness conditions Q, (16), (29), and (30). In our computation we use the software BAMG [9] to 
generate adaptive meshes as M-uniform meshes for given metric tensors. 

Example 5.1 (Positive definite Hessian, mild anisotropy). This example is in the form of BVP gj 
with / and g chosen such that the exact solution u is given by a simple quadratic function 



u(x,y) 



2 i 2 

x + ay , 



a > 0. 



Two different values of a are used, a = 1 and a = 25. 



We first look into mesh conditions ([9j) and (10). Denote the left-hand side of the inequalities by 
Q eq (x) and Qaii(x). A typical palette plot for these functions is shown in Fig. [2j We can see that 
in both cases, 

0.9 < Q eq (x) < 2.1, 1 < Q aH {x) < 1.6, 

which gives C eq = 2.1 and C a u = 1.6. In fact, we found that these bounds also hold for all 
computation with this example, indicating that BAMG [9] does a good job in generating M-uniform 
meshes for a given metric tensor. 

Next we consider the convergence of the finite element solution and the recovered Hessian. Fig. [T] 
shows that the H 1 semi-norm of the finite element error converges at a fir st order rate in terms of 

On the hand, the figure 



2.1 



N~2 as N increases, confirming the theoretical prediction in Theorem 
also shows that the recovered Hessian does not converges to the exact Hessian. 

We now look into closeness conditions Q and ( 16 ). The results for a = 1 and a = 25 are shown in 
the left and right columns of Fig. [3j respectively. In Figs. 3a and 3b e in Q is plotted for adaptive 
meshes obtained with all four Hessian recovery methods. We see that s ~ 0.1 and therefor e Q 
holds for all but DLF method for the case a = 1 when the mesh is almost uniform (cf. Fig. 2a). 
However, the situation is different for the case a = 25 when the solution and therefore the adaptive 
meshes are mildly more anisotropic. In this case, e is not small (e > 1) and therefore Q is not 
satisfied for all meshes obtained with QLS and DLF and some meshes obtained with LLS and WF. 



On the other hand, Figs. 3c and 3d show that the ratio Cr+/Cr- associated with (16), although 
slightly rising with the increasing number of elements, remains 0(1) for a = 1 and a = 25 and 



all four Hessian recovery techniques. Thus, closeness condition ( 16 ) is satisfied by the recovered 
Hessian obtained by these techniques. 

Example 5.2 (General Hessian and strong anisotropy). This example is in the form of BVP 
with / and g chosen such that the exact solution is given by 

u(x,y) = tanh(60x) — tanh(60(x — y) — 30). 

This solution exhibits a very strong anisotropic behavior and describes the interaction between a 
boundary layer along the x-axis and a steep shock wave along the line y = x — 1/2. Moreover, the 
exact Hessian is not uniformly positive definite and regularization is needed in the definition of the 
metric tensor. 

A typical palette plot for the left-hand side functions of ^ and (10) is shown in Fig. 4a 
demonstrating that conditions (Tol) and (10) hold with relatively small C eq and C a u. 
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Numerical results for s and Cr^/Cr- are shown in Fig. |4j We see that, unlike the previous 
example, £>1 and thus, the condition Q which requires e to be small is violated for all meshes in 
the considered range of TV" for all four recovery techniques. Notice that e decreases with increasing N 
but is still O(10 2 ) for the very fine meshes. On the other hand, compared to the previous example, 
the ratio Cr+/Cr- has a larger value (~ 10 4 ) for the current example but, nevertheless, it stays 
bounded with increasing N, confirming that ( |29[ ) and (30) are satisfied by the recovered Hessian. 

The fact that the ratio C r + / Cr- has different values in this and previous examples indicates 
that the accuracy or closeness of the four Hessian recovery techniques depends on the behavior and 
especially the anisotropy of the solution. Fortunately, as shown by Theorems 2.1 and 3.1, the finite 



element error is insensitive to the closeness of the recovered Hessian. The error is shown in Fig. 4e 
as a function of N. 

Finally, Fig. 4d shows that ah, computed through (34), is close to the exact value a defined in 
& 



6 Conclusion and further comments 



In the previous sections we have investigated how a nonconvergent recovered Hessian works in 



mesh adaptation. Our main results are Theorems 2.1 and 3.1 where an error bound for the linear 
finite element solution of BVP §5§ is given for quasi- M-uniform meshes corresponding to a metric 
depending on a recovered Hessian. As conventional error estimates for the H 1 semi-norm of the 
error in linear finite element approximations, our error bound is of first order in terms of the average 
element diameter, N~d l where TV" is the number of elements and d is the dimension of the physical 
domain. This error bound is valid under the closeness condition (16) (or (29) and (|30|)), which is 
weaker than Q used by "Adaptive generation of quasi-optimal tetrahedral meshes" II] and "An 
adaptive algorithm for quasioptimal mesh generation" jl8]. Numerical results in Sect. [5] show that 
the new closeness condition is satisfied by the recovered Hessian obtained with commonly used 
Hessian recovery algorithms. The error bound also shows that the finite element error changes 
gradually with the closeness of the recovered Hessian to the exact one. These results provide an 
explanation on how a nonconvergent recovered Hessian works in mesh adaptation. 



In this work the closeness condition (16) (or (29) and (30)) has been verified only numerically. 



Developing a theoretical proof of the condition for some Hessian recovery techniques is an interesting 
topic for further investigations. 
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Figure 3: Example 5.1 u(x, y) = x 2 + ay 2 . 
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Figure 4: Example 5.2 u(x, y) = tanh(60x) — tanh (60(x — y) — 30). 
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